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Abstract: In this paper, we study the success rate of the reconstruction of 
objects of finite extent given the magnitude of its Fourier transform and its 
geometrical shape. We demonstrate that the commonly used combination 
of the hybrid input output and error reduction algorithm is significantly 
outperformed by an extension of this algorithm based on randomized 
overrelaxation. In most cases, this extension tremendously enhances the 
success rate of reconstructions for a fixed number of iterations as compared 
to reconstructions solely based on the traditional algorithm. The good 
scaling properties in terms of computational time and memory requirements 
of the original algorithm are not influenced by this extension. 
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Inverse problems play an important role in many areas of physics. In an inverse problem we 
intend to reconstruct physical properties ("scattering potential") of some object given the result 
of the interaction ("scattering") with a well-defined input signal. Typically, propagating wave 
fields (e.g., acoustic waves, electromagnetic waves, . . . ) are used to scan the object. However, 
for rapidly oscillating electromagnetic wave fields like x-rays, detection devices are only capa- 
ble of recording the average intensity of the scattered wave field. Its phase information, which 
also contains important information about the scattering potential, is lost during the measure- 
ment process. Despite the lost phase information, x-rays are an important, non-destructive tool 
for the investigation of structural properties of matter because they exhibit high penetration 
depth and high spatial resolution 0J[3). Note, that 3 rd generation synchrotrons [;4] provide bril- 
liant x-ray beams with a very high degree of coherence and with very high flux, so that even 
dot- and wire-like nanostructures can be investigated [2. 5T- [TT1 . 

In this paper, we introduce a new algorithm - the (HIO+OR)+ER-algorithm - which is ca- 
pable of reconstructing the missing phase information in a variety of cases for which the tradi- 
tionally and commonly used HIO+ER-algorithm 11211131 (reviewed in sec. [2]) fails. Before we 
introduce this new algorithm in sec. [2] we shortly review the phase retrieval problem in sec. [1] 
In particular, we want to stress the fact that the (HIO+OR)+ER-algorithm is based on random- 
ized overrelaxation 01411151 and, thus, does not require additional mathematical constraints |fl6l 
compared to the HIO+ER-algorithm. We reformulate our approach using projection polyno- 
mials which facilitate a straight forward investigation of randomization of the HlO-algorithm 
without correlations in the coefficients enforced by overrelaxation. Section [3] discusses numer- 
ical aspects how to measure the convergence of phase retrieval. Finally, sec. [4] illustrates the 
advantages of the (HIO+OR)+ER-algorithm with respect to the traditional HIO+ER-algorithms 



in detail. For this purpose, we focus on three numerical examples each with different spe- 
cific properties and compare the success rate of reconstructions based on our (HIO+Or)+ER- 
algorithm and on the HIO+ER-algorithm. Moreover, we investigate the sensitivity of our al- 
gorithm to the choice of its parameters. In this paper, we focus on the algorithm itself and its 
improved convergence properties for specific, well-known problems related to phase retrieval 
based on the HIO+ER-algorithm. A separate publication will illustrate the improved features 
of the (HIO+OiO+ER-algorithm for the reconstruction of the displacement field in inhomoge- 
neously strained nanocrystals. 

1. The phase retrieval problem 

Consider some unknown object /i n (jc) in ^-dimensional direct space ("position space"). We 
assume that (i) the shape § of f m (x) (i.e., the smallest domain in direct space for which /j n (3c) 7^ 
0) and (ii) the quantity A-? = |FT{/; n (x)}| (i.e., the amplitudes of its Fourier transform) are 
known. It has been proven that this information is sufficient to reconstruct the object /; n (x) if 
the shape § infinite and the dimension d of direct space is at least two l|2l [T7l - [T9l . Physically, 
the second condition is fulfilled in scattering experiments which can be described in lowest 
order Born-approximation (also called kinematical approximation). This is typically true for 
nanostructures illuminated by coherent x-ray beams. 

Suppose we discretize a rectangular region D C W 1 in direct space, which fully contains our 
object /i n (x) (i.e., § C D), with Nj, i 6 {1,2,..., d}, equidistant points. The discrete Fourier 
transform (DFT) maps this set of points to a mesh M in reciprocal space which is related to the 
Fourier transform of our continuous direct space object fi n (x) ll20l . Hence, we have 



where x = {x\ ,X2, ■ ■ ■ ,Xci) T corresponds to the equidistantly sampled points in direct space de- 
fined by D and <t>£ are the phases associated with the amplitudes Aj in reciprocal space. 

In direct space, every point x inside the object domain § corresponds to (up to) two real 
unknowns (magnitude and phase). However, outside the domain S, fm(x) is precisely zero for 
all positions x. Hence, each point outside § defines two real equations for the determination of 
the N real phases <t>j, if we require fm{x) = V x £ S. Therefore, we can hope to reconstruct 
the missing phases <t>£ on the sampled grid M once 



where N§ is the number of points inside the region S. The fraction a = N/N§ is called oversam- 
pling ratio. From our elementary discussion, we can conclude that a — 2 is a lower bound for a 
successful reconstruction, if no additional a priori knowledge is available. For more details on 
oversampling, we refer the reader to the investigations of Elser et Milliane |j2~T| and of Miao et 
al. ll22l (and the references therein). 

With that knowledge, we restate the phase retrieval problem more precisely: Given the shape 
S of some object f m (x) and the modulus A^ of the Fourier transform of f m (x) on a mesh M that 
satisfies Eq. we try to reconstruct the phases <J>j and, thus, the object f m (x). 

Note, that the reconstruction is only unique up to the unavoidable inherent symmetries of 
the Fourier transform |]2][T7][I1], i.e., a global phase shift (typically irrelevant), a plane wave 
modulation e lkx in reciprocal space (fixed by the position of the shape § inside D) and the 
degeneracy between the object fy and its complex conjugated, inversion symmetric object fS. 
This last ambiguity constitutes a severe problem if the shape § is inversion symmetric, but the 




(1) 



2(N-N S )>N <-> N>2N S , 



(2) 



object /;„ does not possess this symmetry 0131231 . For the remainder of this paper, we restrict to 
domains § which are not inversion symmetric. However, even if we exclude objects with inver- 
sion symmetric shape, it is highly non-trivial to formulate a phase retrieval algorithm that works 
automatically without manual supervision and tuning of parameters during the reconstruction. 
The (HIO+OR)+ER-algorithm which we propose in the next section is capable of performing 
this task for a large class of objects f m . 

2. Reconstruction algorithms 

In this section, we will shortly review the traditional combination of the hybrid input output 
and error reduction (HIO+ER) algorithm 111 2111 31 and propose our extension to this algorithm. 
Throughout the entire reconstruction process, we assume that we know the exact geometrical 
shape § of the object which we reconstruct. 

2.1. Hybrid input output and error reduction 

A fast and efficient iterative phase retrieval algorithm is based on alternating projections of 
a trial solution onto the constraints in direct space and in reciprocal space. Hence, a single 
iterative step (i) of this algorithm, which is called error reduction (ER) algorithm [ 1 2\ , is defined 
by 

ft 1] =PsPAf^H BR f^ . (3) 

Typically, initial phases 4>i° are chosen randomly. Here, P§ and Pa are projection operators in 
direct space and reciprocal space respectively, i.e., 

and 

P^=A|«p(iaig(^°j) . (4b) 

For notational simplicity, we assume that any operand of Pa or P§ is transformed to the proper 
space (i.e., Fourier transform to reciprocal or direct space respectively) before the operator Pa 
or Pg is applied (e.g., Pa/% means Pa FT{/i'' ) }). 

The projection operator Pa is non-linear, non-convex and non-unique 11241 . The values \g~ | 
are important for determining the difference to the (experimentally accessible) input data Aj 
(see Eq. (fT8l)). The ER-algorithm is a local minimizer of a suitable chosen error metrics H12II131 . 
However, in practice, phase retrieval problems involve many local minima. Hence, the ER- 
algorithm will in general not converge to the correct solution f m (x), but to some local minimum. 
In addition, the error metric may stagnate for many iterations before decreasing further. More 
information on stagnation problems and the ER-algorithm can be found in 1 1 2 131 1231 . 

Due to these problems, further algorithms, which try to avoid stagnation and convergence 
to local minima, have been developed. A very important algorithm is the hybrid input output 
(HIO) algorithm proposed by Fienup 112111311231 . This algorithm is also an iterative procedure 
which is defined by the map 

where j5 is a real parameter typically chosen in the range [0.5; 1.0] 1231 . Note, that this definition 
can be rewritten as 

/( ,+1) = [l-P s -)3PA + (l+ i 8)PsPA]/i , ' ) ^i/Hio(j3)/i ,) . (6) 



The convergence properties of the HlO-algorithm have been investigated in more detail in 
C2|25]|26). 

This algorithm is much stronger in avoiding stagnation and local minima. However, the HIO- 
algorithm shows its full potential only if it is combined with the ER-algorithm. One step of this 
combined HIO+ER-algorithm consists of Nhio iterations of the HlO-algorithm followed by 
Ner iterations of the ER-algorithm. This combination is more successful than both algorithms 
on their own as it was observed in practice lfT3ll23l . Although this algorithm is already very 
powerful, it is not yet satisfactory for many objects as we will show in sec.H] 

2.2. Hybrid input output with randomized overrelaxation ((HIO+Or)+ER) 

The HIO+ER-algorithm can be further improved by including randomized overrelaxation lfl4l 
15 27 ). Basically, overrelaxation corresponds to replacing a projection operator by 

0^ = 1 + ^^-1), (7) 

where A^ is a real constant called relaxation parameter lfT4l . The limiting case Q^-x^ — > P/j is 
obtained for A,^, — ^ 1 . 

Overrelaxation (without any randomization) has been investigated for convex problems lfT5l 
and in connection with the ER-algorithm for phase retrieval lfl4l . Moreover, overrelaxation is 
also included in the difference map algorithm proposed by Elser in 11271 . In the difference map 
algorithm with overrelaxation, the iterative step for our set of constraints is given by l28l 

= [1+J3 (PsQ A ;X A -PaQs-m)] fP = /W/U A ,A s )/f , (8) 

where Elser proposes to choose the parameters as Aa — A§ = J3 _1 11271 . A very nice comparison 
of several phase retrieval algorithms and a benchmark thereof can be found in l28l . Marchesini 
demonstrated in [29 1 that for his particular numerical example an additional low-dimensional 
subspace saddle-point optimization was also able to overcome stagnation of the traditional 
HIO+ER-algorithm. In our extension, such an additional optimization is not necessary. 

We propose to replace the (non-linear and non-convex) projection operator Pa in reciprocal 
space in the HlO-algorithm itself by its relaxed expression 

Qa ; a a = 1 + ^(Pa-1) • (9) 

The direct space assembly in Eq. (O of the HlO-algorithm remains unchanged. This corre- 
sponds to changing Eq. (O to 

ft l) = [1 - Ps- PQam + (1 +P)PsQaaJ $ = Hmo+o R {PM)f! ] ■ ( 10 ) 
If we rewrite this expression in the projectors P§ and Pa, we get 

H m o+o R (0 , Xa) = [1 + J3 (X A - 1)] + [j3 - X A - pX A ] P§ - pX A P A + [(1 + j3)A A ] P§P A .(11) 

The deviation /3 (Xa — 1 ) from the identity operator in the first term can neither be represented 
by the traditional HlO-algorithm for any value of j3 (see Eq. ©) nor by the difference map 
algorithm for any (j3, Xa, A§) (see Eq. dH)). In both cases, the previous iterative result /- is 
weighted with 1 or projected at least once in either direct (P§) or reciprocal space (P^) before 
being included in the calculation for/i' +1 '. 

Finally, we have to choose suitable values for the additional parameter Aa in Eq. (0. If we 
restrict to a fixed set of values Aa for all iterations, we risk simply exchanging one trap by 
another or one local minimum by another. Trying to overcome local minima and stagnation, 
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Fig. 1. Graphical illustration of the (HIO+Or )+ER-algorithm and its building blocks. 



we propose a new approach: For each iteration, we randomly select the parameter Xa- More 
precisely, for the remainder of this paper, we investigate a uniform random distribution in the 
range [1 — V, 1 + v], V > 0, for Xa which is reassigned each iteration. Unless stated otherwise, 
we choose v = 0.5. Note, that a fixed value Xa = 1 for all iterations corresponds to the usual 
HlO-algorithm. In order to distinguish our new extension from the traditional HlO-algorithm, 
we call our extension HIO+OR-algorithm. The power of this extension is strongly supported by 
the fact, that the HIO+OR-algorithm is capable of reconstructing objects without including the 
ER-algorithm with significant success rates. If the HIO+OR-algorithm is combined with ER, 
we call the algorithm (HIO+OR)+ER-algorithm (see Fig. |TJ. 

If we look at Eq. (fTTT i from a bird's eye view, the most important difference of the operator 
Hhio+o r (P,Xa) to the other algorithms which have been proposed {28] is the fact, that the 
traps and tunnels lfl4l of the operator in general depend on Xa and, therefore, should change 
from iteration to iteration as a consequence of randomization. Thus, stagnation in a specific trap 
or tunnel is strongly reduced, as most traps and tunnels are no longer persistent throughout the 
iterative reconstruction process. However, the true solution f m (x) is a fixed point of the iterative 
procedure defined in Eq. (fTTT i for all values of Xa- 

2.3. Projection polynomials 

The extension of the HIO+ER-algorithm to the (HIO+OR)+ER-algorithm is based on two in- 
gredients: overrelaxation and randomization. In order to gain a deeper understanding of the 
(HIO+OR)+ER-algorithm, it is useful to investigate both ingredients separately. Studying the 
HlO-algorithm extended by overrelaxation of P^, but without randomization is straight for- 
ward: we simply keep a fixed, predefined overrelaxation parameter Xa throughout the entire 
iterative reconstruction process. In this section, we present a framework for randomization of 
iterative projection algorithms performed in a way which is not based on overrelaxation of P^. 
Although this framework can be applied for studying randomization of any iterative projec- 
tion algorithm, we limit our numerical simulations (in section [U to parameter values whose 
deterministic contribution coincides the traditional HlO-algorithm. 



For this purpose, consider the projection polynomial operator 



#ProjO,4,4)=M+£ [cs, 2 „(PsPa)"+Ca,2„(PaPs)"] 



11=1 



odd 
"Max 



E [c S ,2, 1+ iP§(PaPs)"+ca,2, 1+ iPa(P§Pa)"] , (12a) 

for obtaining the next iterative solution in the manner of Eq. © or Eq. (fTTT i. «^a" * s gi ven by 
the largest integer smaller or equal to | . «^ax ^ s determined by the largest integer smaller or 
equal to p is the maximum number of successive projection operators which should be 
included in Hp mi . 

This operator exploits the fact, that any combination of projection operators P'^ j . . . P^'" of 
two different kinds <fjj £ j £ {l,...,m}, reduces to one of the four building blocks 

Ps (PaP§)"» Pa (P§Pa)", (P§Pa)" and (P A Ps)" for some integer n > 0. This is a conse- 
quence of the defining property P| = P^ (idempotence) of projection operators (which implies 
P| = P| V n > 1). Note, that the product of two non-commutative idempotent operators (like 
Pa and P§) is no longer idempotent. The special case of the identity operator 1 and a single pro- 
jection operator is included for the case « = in those building blocks. However, the identity 
operator is included twice, namely for n = for (P§Pa)" and (P^P§)". Therefore, one spuri- 
ous coefficient in the projection polynomial is introduced if both building blocks are included 
for « = 0. This is avoided by restricting those two building blocks to n > 1 and including the 
identity operator separately. Therefore, Eq. (I12at represents the most general polynomial ex- 
pression which contains a maximum of q successive projections of two different kinds, if we 
incorporate the idempotence of projection operators. 

In order to guarantee, that the true solution f m (x) is a fixed point of Hp m] for any choice of its 
parameters (£>,c§,ca), it is necessary to enforce one constraint on the coefficients: If the input 
for the projection operators P§ and Pa coincides with the true solution f m (x), the projection 
operators P§ and Pa reduce to the identity operator 1 . Therefore, given the true solution /; n (x) 
as input, the operator //p ro j simplifies to the identity operator 1 for any choice of its parameters 
(^,c§,ca) if and only if the additional constraint 

p 

fo = l-£ [c„ s + c„a] ■ (12b) 
n=i 

is fulfilled. Consequently, for a maximum of p successive projections in the projection poly- 
nomial Hp ro j, only 2p free parameters appear in H Pro j by incorporating the idempotence of 
projection operators. 

Given this fix point property for the true solution fm(x), we can perform proper randomiza- 
tion of this operator. For this purpose, we split our coefficients c^ n , £, £ {S,A}, n > 1, in a 

deterministic part c^°j and a random part r^c^ and assume r^ n to be uniformly distributed 

in the range [—1,1], i.e., 

c 5,"= c S +r S>" C 8 ■ (13) 

(R\ 

If we set some coefficients n ^ and draw statistically independent random values for r^ n , 
we can investigate randomization of the HlO-algorithm without the correlations in the coeffi- 
cients „ introduced by overrelaxation. b is fixed by Eq. (1 12bb and, thus, is fully correlated to 
the values of n . 



In order to obtain the traditional HlO-algorithm in the limit n — > for all £, and n, we have 
to choose (at least) q = 2 (see Eq. ©). Therefore, flproj reduces to 

Hmo R {b, c§a , cs,i, ca,i, c A ,2) = b 1 + c §il P§ + ca.iPa + c§, 2 P§P/i + c a ,2PaP§ (14) 

with the additional constraint b = 1 — E«=i [ c n.S + c„A]- The deterministic contribution c^°j 
reproduces the traditional HlO-algorithm (see Eq. (|6j) if we set 

eg? 1, cg = l+j8, cg=0. (15) 

After replacing P^ by Q A ;X A ( see Eq- ©X the coefficients cp „ are parametrized by two param- 
eters j3 and X A as 

c s,i =j3 -X A -fSX A , c a ,i = -/Ua, cs,2 = (1 + J3)Aa , c A ,2 = 0, (16a) 
which corresponds to 

c s ,i = -1-7a(1+J3), ca,i = -jS(1 + 7a), c Sj2 = (1 + /J)(1 + y A ) , c A , 2 = 0, (16b) 

if we substitute X A = 1 + y A . y A is uniformly distributed in [— V,v]. The deterministic contri- 
bution in Eq. ( 116bt is identical to Eq. ( fT31 ) and reproduces the traditional HlO-algorithm in the 
limit v —> 0. However, the random contribution of each coefficient n is determined solely by 
the value of y A and, thus, not statistically independent, but fully correlated. This is in strong 
contrast to a statistically independent randomization of each coefficient c^ n separately. 

Consequently, we can apply the framework of projection polynomials to investigate whether 
the benefits of randomization rely on the precise correlations evident in Eq. ( |16b| ) or if random- 
ization of the traditional HlO-algorithm implies significant benefits also in absence of these 
correlations. 



3. Monitoring convergence of the reconstruction 

In order to investigate the success and convergence properties of the reconstruction, we monitor 
three parameters. First of all, we monitor the change of the reconstructed object /w from itera- 
tion (z — 1) to iteration (/). Second, we measure the mathematical distance to the (sampled) true 
solution fm(x). And last, but not least, we measure the deviation in reciprocal space between 
the a priori given reciprocal input data A- k and the magnitude of the Fourier transform of the 
reconstructed object. All three parameters yield different information: 

The distance of the reconstructed object /W after iteration (/) to the true direct space object 
/in is, of course, the best measure for the quality of the reconstructed image. We define this 
mathematical distance by an angle 



= arccos 



(| (/ (0 ; /in )|V</« ;/('■)> </ in ; ~f~) \ . (17) 



The absolute value in the nominator of the argument of the arccos eliminates the influence of 
the undetermined global phase in /('). Moreover, <pW is identical in direct and reciprocal space 
(invariance of scalar product upon unitary basis transformations like the Fourier transform). 

A dimensionless, normalized error measure eW which depends only on the a priori known 
(measured) amplitudes and not on the true direct space solution / n (x) is 

L - — %\$H*)\-m) , (is) 



(A; A) \\A\ 



2 




(a) This object belongs to Eq. \22\ . (b) This object belongs to Eq. (23). 



Fig. 2. Pure phase objects used for the investigation of the convergence of the 
(HIO+OR)+ER-algorithm. The magnitude of both objects is constant. The phase is plot- 
ted using a HSV color-bar, however, the region outside the shape § has been set to black. 
The oversampling ratio is a = 8.456. 



where ||-|| is the p-norm. Since the construction of the absolute value is a non-linear, non- 
invertible map, the result can in general be quite different from the angle defined above. In fact, 
two quite different direct space objects f m (hence quite different complex Fourier transforms) 
may posses an extremely similar distribution of the magnitude A in reciprocal space H16II181 . 

In addition, from a numerical point of view, the change of the reconstructed object /W from 
iteration to iteration is very important, too. This must not be confused with the change of the 
non-invertible map eW from iteration to iteration: Moving along a connected line of constant 
error metric does not change the error metric itself, but the reconstructed object /W may change 
arbitrarily. Nevertheless, the algorithm may converge (i.e. no change from iteration to iteration), 
but not to the true solution /;„ of the problem II14II16I . For simulated data, we can judge if the 
algorithm converged to the true solution f m by Eq. dPTl . A suitable choice for observing the 
change from iteration to iteration is given by the angle 
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In contrast, the p-norm 5™ = /( ! ~ _1 ) — fw of the difference is not optimal for measuring 

p 

the convergence of the reconstruction process, because it does not eliminate global phase shifts 
from iteration (i — 1) to iteration (;). 

4. Numerical examples 

In this section, we demonstrate the improvements resulting from randomized overrelaxation for 
three carefully chosen numerical examples: First, we consider two pure phase objects 

/(x)=exp(i^(x))n s (x) , (20) 

where £ (x) is a real function and £2g(x) is the shape function 

_ _ J ifx£§, .... 
%(*) = | ! iffe g (21) 

The reconstruction of the first phase object (see Fig. |2(a)) is plagued by (at least) one strong 
local minimum far away from the true solution f m which results in severe problems for the 





Fig. 3. Investigation of the (HIO+OR)+ER-algorifhm for a purely real object f m with strong 
variation of its magnitude over short length scales (Oversampling ratio a = 5.06). In Fig. 
|(b)[ continuous lines represent the (HIO+OR)+ER-algorithm, isolated dots the HIO+ER- 
algorithm. A pure HIO+OR-calculation without ER is included as black, dash-dotted curve. 



traditional HIO+ER-algorithm. The reconstruction of the second phase object (see Fig. |2(b)| i 
is plagued by many local minima near the perfect solution /;„, in particular, by stripes ll23l . 
Both phase objects have a smooth variation of their phase, i.e., phase variations of 2n extend 
over several sampling points. As a third test object, we investigate a purely real object f m {x) 
with strong variation of its magnitude over very short length scales and a different shape (see 
Fig. |3(a)| i. In fact, the success rate of the traditional HlO-algorithm without randomized over- 
relaxation was lowest for this third example. From the numerical results, we conclude that 
randomized overrelaxation is a powerful extension for a wide class of objects and not only for 
objects possessing some particular features. 

The first pure phase objects which we used in our numerical simulations is given by 



(22) 



where £,\ is defined according to Eq. (l2Qb and the parameters were chosen as b\ = 1.5515 and 
c\ = —1.835. The shape § has been restricted to the triangular region defined by (0,0), (0,a) 
and (a,0) with a = 0.97. The object has been sampled on an equidistant grid with Ni x N2 = 
256x256 data points in the interval (x,y) = (—1,-1) to (x,y) = (1,1). The resulting object is 
depicted in Fig. |2(a)| 

The second object is defined by 



&(jc,y) = (2n) 



{x-b 2 f + (J - C 2 ) 2 + ^273 



(23) 



with the same shape § and sampling parameters as in the first example. The parameters were 
set to £>2 = —0.3515 and C2 = 0.535. The phase field in this case is depicted in Fig. |2(b)| 

A reconstruction is considered successful once the angle to the input object (see Eq. dTD i) 
falls below some given limit ^"converged (=0.05° unless stated otherwise). We repeat the recon- 
struction process for Nn e si initial trials (random phases at each fc-point, amplitudes A^ equal 
to given amplitudes At). Finally, we calculate the success rate which tells us the percentage of 
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(a) First phase object denned in Eq. i22\ . 
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(b) Second phase object defined in Eq. j23\ . 



Fig. 4. Comparison of the success rate of reconstructions of pure phase objects (see Eq. 
( f20t . l l22l > and l[23j) with the HIO+ER- and the (HIO+0 R )+ER-algorithm. The parameter j8 
was fixed to 0.85. Continuous lines represent results of the (HIO+OR)+ER-algorithm, iso- 
lated dots of the HIO+ER-algorithm. A pure HIO+OR-calculation without ER is included 
as black, dash-dotted curve. 
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Fig. 5. Long-term stagnation of the success rate of the traditional HlO-algorithm (with- 
out overrelaxation and without randomization) for the first phase object (see Eq. 122b for 
different choices of the internal parameters (A'hiOi^'er)' but fixed = 0.85. 



reconstruction processes that have been successful up to iteration (/). A good reconstruction 
algorithm 

(i) should reach a success rate of almost 100%, 

(ii) should not depend on the starting point (as long as no good starting point is available), 

(iii) should perform the reconstructions with little computational effort and 

(iv) should possess these properties for a wide range of its internal parameters (i.e., j3, V, A^hio 
and Neb, in case of the (HIO+OR)+ER-algorithm). 

Each iteration in the HIO+ER-algorithm and its (HIO+OR)+ER-extension can be imple- 
mented very efficiently exploiting the good scaling properties of the FFT-algorithm. Therefore, 




Fig. 6. Investigation of the sensitivity of the (HIO+OR)+ER-algorithm with A'hio = 50 and 
A'er = 20 on the choice of the parameter p and on the range of the uniform distribution 
determining the relaxation parameter A4 for the first phase object (see Eg, 122), 
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(a) First phase object defined in Eq. [H\ . 




Iteration [ J 



* X =0.91 


O X A =0.99 


x x A =i.oi 


X A =0.9iS 


X =1.00 

A 





(b) Second phase object denned in Eq. )23t . 



Fig. 7. Success rate of the HIO+ER-algorithm for fixed overrelaxation (no randomiza- 
tion). Parameters are chosen as (A'hio = 130JVer = 10) and j8 = 0.85. 



as long as the number of iterations is not too high, the third requirement is clearly met by both, 
the HIO+ER and (HIO+0 R )+ER-algorithm. 

In order to investigate the other points, we plotted the respective success rates for the 
HIO+ER- and (HIO+OR)+ER-algorithm for different internal parameters in Fig. 5] We nor- 
malized the computational effort, i.e., one iteration of the (HIO+OR)+ER-algorithm with 
Nmo = 130 and A%r =10 corresponds to two successive iterations with A^hio = 50 and 
Ner = 20. This prefactor of 2 can be found in the legend of the graph in front of the square 
brackets for each calculation. 

The inclusion of overrelaxation improves the success rate for the first phase object defined 
by Eq. (l22l) for each combination of A^hio and A%r: in contrast to the HIO based calculations, 
it quickly converged to 100%. In fact, overrelaxation is so powerful that stagnation in the HIO- 
algorithm is eliminated entirely and intermediate ER-iterations are no longer necessary (see 
black dash-dotted curve in Fig.|4|i. The success rate of the traditional HIO+ER-algorithm stag- 
nates on some level far below 100%. This stagnation is significant in the long term behavior, 



as we illustrate in Fig. E]by plotting the success rate of the traditional HIO+ER-algorithm for 
50 times as many iterations as in Fig. |4(a)| Note, that all trials that did not converge after 2500 
iterative steps using the traditional HIO+ER-algorithm with Nmo = 130 and A^r = 10 in Fig. 
|5]are separated from the true solution by an angle greater than <p = 55°. 

If we consider the second phase object defined by Eq. ( 1231 . the traditional HIO+ER- 
algorithm performs better here than for the first phase object. For the one specific set of pa- 
rameters A^hio = 130 and Ner = 10, the success rate of the traditional HIO+ER-algorithm is 
even slightly superior to the (HIO+OiO+ER-extension, but for this specific set of parameters, 
the success rate converges to 100% very quickly for both algorithms. However, for the set 
of parameters (Nhio = 50, Ner = 20) and (Nhio = 40, Ner = 30), the success rate quickly 
converges to 100% only for the (HIO+OR)+ER-algorithm, but not for the HIO+ER-algorithm. 
Hence, the (HIO+OR)+ER-algorithm is much more stable with respect to the particular choice 
of the internal parameters A^io an d Ner including, in particular, the case Ner = 0. 

The traditional HIO+ER-algorithm also has severe problems performing a successful recon- 
struction for the real object depicted in Fig. |3(a)| The success rate for the reconstruction process 
is very poor if we apply our criterion Pconvereed < 0.05° to determine successful reconstructions. 
The (HIO+OR)+ER-algorithm is much more successful and reaches success rates of 100% with 
few iterations (see Fig. |3(b)| i. Furthermore, reaching the success rate of 100% does not sensi- 
tively depend on the internal parameters A^hio and Ner- Again, keeping Ner small compared 
to A^hio or eliminating ER-iterations completely yields the best results. Keep in mind, that nei- 
ther the reality constraint nor the positivity constraint of the object have been exploited during 
reconstruction. 

For the second and third test object, all trials that failed to converge to the true solution f m 
basically evolved towards the correct direction. However, at some point, stagnation in form 
of stripes close to the true solution f m appeared. These stripes are very persistent and prevent 
the usual HIO+ER-algorithm from further convergence to the true solution /;„. Typically, the 
mathematical distance <p of the objects containing stripes to the true solution f m is in the order 
of 0.05 to 5.0 degrees. In case of the purely real test object, after 150 iterations all reconstructed 
objects reached a distance (p < 0.5° (= 10 • ^converged) for HIObo — ERjo and 3 1 % for HIO50 — 
ER?o- The remaining 71% of the objects for HIO50 — ER20 reached a distance <p < 5° (= 100 ■ 
(^Converged) using traditional HIO+ER-algorithm and 150 iterations. Note, that the change #W 
(see Eq. ( fT9l )) from iteration to iteration does not vanish, i.e. numerically, the algorithm did not 
converge within the given number of iterations. 

Figure [6(a)| shows the behavior of the success rate of the (HIO+OR)+ER-algorithm for differ- 
ent values of the parameter v. In addition, Fig. |6(b)| depicts the success rate for different values 
of the parameter j3. Within the range j3 € [0.5, 1.0], which is typically used in the traditional 
HIO+ER-algorithm, the (HIO+OR)+ER-extension does not sensitively depend on the particu- 
lar choice of /3 . Therefore, the (HIO+OR)+ER-algorithm does not depend sensitively on any of 
its internal parameters A^rio, Ner, P and v. Hence, it fulfills all four requirements on a good 
reconstruction algorithm which we stated above. 

Two concepts (overrelaxation and randomization) have been exploited to modify the HIO- 
algorithm: Fig. |7]depicts the success rate, if overrelaxation is performed with a static parameter 
Xa- We see extremely sensitive behavior of the success rate on the precise value of Xa if it 
is kept fixed. Whereas increasing up to 1.02 improved the success rate for the first object, 
it completely vanished for Xa = 1-01 for the second phase object. For a deviation from 1.00 
greater than ten percent, almost no successful reconstructions have been observed for any of 
the three test objects. In addition, we do not a priori know the (non-universal) optimum value 
for Xa which depends on the object f m . Consequently, static overrelaxation does not fulfill the 
requirements for a good reconstruction algorithm stated above. 
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Fig. 8. Comparison of the success rate of both frameworks which provide a generalization 
of the traditional HlO-algorithm based on randomization. Continuous lines illustrate the 
behavior for randomized overrelaxation of (see Eq. l|10|>). whereas dots represent the 
behavior of the success rate resulting from independent randomization of the coefficients 
in a projection polynomial (see Eq. d!3t and Eq. d!4t). In both cases, the deterministic 
contribution is equivalent to the traditional HlO-algorithm with parameters Nmo = 140 
and p = 0.85. No ER has been performed. 

Based on the discussion in sec. 12. 31 we test randomization of the HlO-algorithm for j3 =0.85 
without employing overrelaxation by setting the amplitudes for the random contribution to 

the coefficients n in the operator H\\\o R defined in Eq. ( [Pfl ) equal to e g j = = Cg^ = 

C A2 = ^ = 0-2. The deterministic contribution is set according to Eq. (TT3T >. The result 
for all three test objects is illustrated in Fig. [8] and compared to the result for the HIO+Or- 
algorithm. Clearly, pure randomization without overrelaxation performs equally well as the 
HIO+OR-algorithm for all three test objects. However, in our opinion, the HIO+Of>-algorithm 
should be preferred over randomization without overrelaxation, because (i) it requires only 
two Fourier transformations for each iteration (lower computational effort) and (ii) it has less 
degrees of freedom (one uniform random distributions instead of four). Note, however, that due 
to the large number of degrees of freedom of projection polynomials, further optimization of 
this approach is possible (including more advanced statistical correlations in the coefficients 
of the projection polynomial). In addition, more elaborate random distributions might improve 
both approaches further. 

In summary, we succeeded in overcoming stagnation of the traditional HIO+ER-algorithm 
for various objects. For this purpose, the concept of overrelaxation was applied to the recip- 
rocal space projection in the HlO-algorithm which produced additional degrees of freedom. 
Randomization of these additional degrees of freedom from iteration to iteration was exploited 
to prevent stagnation in local minima, especially in traps or tunnels. In particular, this was 
demonstrated for two pure phase objects with smooth phase variation (with different behavior 
with respect to the traditional HIO+ER-algorithm) and a purely real object with strong am- 
plitude variations over short length scales. The improved algorithm is far less sensitive to the 
specific choice of the internal parameters Nmo an d -^er than the traditional HIO+ER-algorithm. 
Moreover, the (HIO+OR)+ER-algorithm remains almost insensitive to the particular value of 
the internal parameter j3 in the range j3 G [0.5, 1.0]. Furthermore, the the good scaling prop- 
erties in terms of computational time and memory consumption of the HIO+ER-algorithm are 
fully maintained in the (HIO+OR)+ER-algorithm. In conclusion, a reconstruction based on the 
(HIO+OR)+ER-algorithm exhibits significantly enhanced stability and success probability in 
comparison with the traditional and commonly used HIO+ER-algorithm. 



